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We showed that a wetting layer in epitaxially strained thin films which decreases with increasing 
lattice mismatch strain arises due to the variation of nonlinear elastic free energy with film thickness. 
We calculated how and at what thickness a flat film becomes unstable to perturbations of varying 
size for films with both isotropic and anisotropic surface tension. We showed that anisotropic sur- 
face tension gives rise to a metastable enlarged wetting layer. The perturbation amplitude needed 
to destabilize this wetting layer decreases with increasing lattice mismatch. We also studied the 
early evolution of epitaxially strained films. We found that film growth is dependent on the mode 
of material deposition. The growth of a perturbation in a flat film is found to obey robust scaling 
relations. These scaling relations differ for isotropic and anisotropic surface tension. 
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I. INTRODUCTION 

The growth of epitaxially strained thin films in which 
there is lattice mismatch between the substrate and the 
film is of major importance in the fabrication of semicon- 
ductor and optoelectronic devices. The lattice mismatch 
generates strain in the deposited film, which can cause 
film instability unfavorable to uniform flat film growth. 
The strained film can relax either by the introduction of 
dislocations or by the formation of coherent (dislocation- 
free) islands on the film surface via surface diffusion. 
These coherent islands can self-organize to create peri- 
odic arrays which can be utilized to create quantum dot 
structures of electronic significance. Understanding and 
predicting strained thin-film evolution is important for 
the improved fabrication of semiconductor devices. 

Early film growth tends to occur via coherent island 
formation as there is an energy barrier to the introduc- 
tion of dislocations. Dislocations occur at island edges 
once islands reach a certain size as the large stress at 
island edges provides a pathway for dislocation forma- 
tion Q. We only consider dislocation- free films since 
many experiments show an absence of dislocations (see 
e.g especially in early film evolution. 

The current work sheds light on the following two prob- 
lems: First, it has been observed experimentally that 
dislocation-free flat films of less than a certain thick- 
ness (the critical wetting layer) are stable to surface per- 
turbations, while thicker films are unstable P-[lC|]. The 
thickness of the wetting layer is substance dependent and 
decreases with increasing lattice mismatch strain [p|dlO|, 
e = (a s — a,f)/a,f, where a s and a/ are the substrate 
and film lattice constants. Above the critical wetting 
layer, 3D coherent islands form. Despite considerable ef- 
forts the physics of the critical wetting layer is poorly 
understood. Namely, why is there a critical, stable wet- 
ting layer and what controls its thickness? As in most 



cases heteroepitaxial growth is done below the roughen- 
ing transition, how does anisotropic surface tension affect 
the thickness of the critical wetting layer? The second 
question we address is how does continuous material de- 
position affect the early evolution of thin films? 

Previous works about the existence and nature of the 
critical wetting layer can be split into two categories: 
those which looked at the dynamic stability of a fiat film 
to small perturbations Jll],[l2| , and those which looked at 
whether a flat film is energetically favourable to a film 
with fully formed faceted islands |l3| - |l6[ . 

The research in the first category addressed substances 
with isotropic surface tension. In these works physical 
parameters (lattice mismatch or surface tension) which 
differ between the substrate and film were smoothly var- 
ied over the substrate-film interface in order to avoid non- 
analyticities. The effect of such smoothing was to create 
a wetting layer. Whilst the choice of a smoothing length 
of the order of a lattice parameter is physically reason- 
able, none of these works gave a physical explanation for 
the smoothing of material parameters over the interface 
or tried to physically calculate the smoothing length or 
form of the transition. 

The research in the second category typically used 
physically motivated methods in order to determine the 
free energy of a flat film of varying depth. Tersoff Ej| 
and Roland and Gilmer fl6|| both used empirical potential 
methods to determine the chemical energy of flat films 
of Ge/Si(001). Tersoff @ then compared this chemi- 
cal potential with that of a bulk strained Ge in order to 
determine whether it is preferable to form islands and 
to predict a wetting layer of 3 monolayers. Roland and 
Gilmer |l6f| saw some evidence of clustering in thicker 
films in molecular dynamics simulations. Daruka and 
Barabasi used an expression for the free energy of a 
flat film of varying depths which fits the results of the 
earlier works |l3|,|ll| in order to compare the energies of 
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flat films and films with fully faceted islands whose en- 
ergies were calculated using continuum elasticity. They 
saw a wetting layer which increased with decreasing lat- 
tice mismatch. Wang et al. jl5| used ab initio methods 
in order to determine the formation energy of flat films 
of varying depths of InAs/GaAs(100). They compared 
this energy with that of a thinner film with fully faceted 
islands whose energies were calculated using continuum 
elasticity. All the above works did not study the issue 
of when a flat film becomes unstable to small monolayer 
perturbations or the dynamics of growth. 

In this paper we show that the variation of nonlinear 
elastic free energy with film thickness can give rise to 
a wetting layer which decreases with increasing lattice 
mismatch strain. We show how and at what depth a flat 
film becomes unstable to perturbations of varying size 
for films with both isotropic and anisotropic surface ten- 
sion. This provides a more realistic estimate of critical 
wetting layer thickness than the studies described above, 
for films in which islands grow from small surface pertur- 
bations rather than being immediately nucleated on the 
flat film. This mode of growth has been seen in many 
experimental systems 0-^,0, especially for films with 
small lattice mismatch, e < 2.5%. As discussed below we 
study the evolution of these small perturbations and ob- 
serve island faceting. We show that anisotropic surface 
tension gives rise to a metastable enlarged wetting layer. 
The perturbation amplitude needed to destabilize this 
wetting layer decreases with increasing lattice mismatch. 

The effects of material deposition on early thin film 
evolution was addressed by Chiu and Gao who looked 
at the evolution of strained films with isotropic surface 
tension when material deposition is constant in the di- 
rection perpendicular to the film surface, corresponding 
to liquid phase epitaxy. In the present paper we look at 
thin film growth when material deposition is constant in 
the direction perpendicular to the film surface and when 
deposition is at a steady rate in the y direction (verti- 
cal to the interface between the film and the substrate), 
corresponding to any directed deposition (e.g, molecu- 
lar beam epitaxy). The latter is a much more common 
method of material deposition in strained film growth. 
We found that the type of evolution seen depends on the 
direction of material deposition. When the deposition is 
constant in the vertical y-direction, the film evolves ac- 
cording to the linear evolution equation, even after the 
surface is no longer a sine function and cusp formation 
occurs. When deposition is constant perpendicular to 
the surface, cusp formation is slowed down at very high 
deposition rates and the surface shows signs of reaching a 
steady-state morphology. We also studied thin film evo- 
lution for faceting films, and found robust scaling laws 
for film growth. 

The rest of the paper is organized as follows. In Sec. 
II we formulate the problem. In Sec. Ill we present the 
general results of linear stability analysis. In Sec. IV we 



describe how we calculated the variation of the nonlinear 
elastic free energy of a flat film with film thickness. It is 
this variation which gives rise to the wetting layer. Cal- 
culations were carried out using a ball and spring model 
in order to determine general qualitative behaviour. Sec. 
V describes the results of the numerical simulations of 
thin film evolution without material deposition, and in 
particular the metastability of the wetting layer. Sec. VI 
describes the results of the numerical simulations of thin 
film growth with material deposition. Some of the re- 
sults on thin film growth without deposition and a brief 
description of the ball-and-spring model for calculating 
the free energy of a flat film appear in our earlier paper 

II- 
ii. PROBLEM FORMULATION 

We model the evolution of a thin film on a substrate 
using continuum theory. The lattice mismatch between 
the film and the substrate creates a strain in the film, 
e . Both the substrate and the film are assumed to be 
elastically isotropic with the same elastic constants. The 
surface of the solid is at y — h(x, t) and the film is in 
the y > region with the film-substrate interface at 
y = 0. The system is modelled to be invariant in the 
z-direction, and all quantities are calculated for a section 
of unit width in that direction. This is consistent with 
plane strain where the solid extends infinitely in the z 
direction and hence all strains in this direction vanish, 
i.e. e xz — e yz — e zz = 0. We assume there is no material 
mixing between the substrate and the film. 

All the results mentioned in this paper relate to vicinal 
surfaces with a very small miscut angle in the z direction 
(see Fig. |l|). Experimentally, surfaces often have such a 
small miscut, as it is very difficult to grow a perfect facet. 
In such surfaces there is no finite energy barrier for the 
formation of an infinitesimal perturbation on the surface, 
since such a perturbation involves only step motion and 
bending, and no nucleation of new steps. For a faceted 
surface with no miscut in the z direction there is a fi- 
nite energy barrier for the formation of an infinitesimal 
perturbation associated with nucleation of step pairs. 

The continuum approximation in the lateral direction 
(x — z plane) is valid as the film is infinite in the z di- 
rection and the smallest lateral surface features we study 
have widths of tens of atoms. However in the vertical y 
direction the films we are studying are sometimes only 
a few monolayers thick. Is the continuum model valid 
for such a film? Can inherently discrete system proper- 
ties, such as the change in free energy as a monolayer is 
added, be interpolated to films of a non-integer number 
of monolayers? 
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FIG. 1. Perturbations of vicinal and fully faceted surfaces. 
The dotted lines represent the cross-sections taken in the xy 
plane which are shown in the bottom graphs. The figure 
shows that nucleation of new steps is needed in order to per- 
turb a facet, but not in order to perturb a surface vicinal in 
the z direction. 

First, consider fully faceted films. As can be seen in 
Fig. [j] there is a qualitative difference between films with 
integer and non-integer numbers of monolayers, and so 
discrete system properties cannot be interpolated to films 
with a non-integer number of monolayers. Hence for such 
films continuum models are not expected to give accurate 
results. To accurately model fully faceted film growth, 
the discrete nature of the steps needs to be taken into ac- 
count using atomistic models and simulations. Modeling 
large systems in such a manner is currently beyond com- 
putational capabilities. However, as can be seen in Fig. 
|l|, the above arguments do not apply to vicinal surfaces. 
In this case a perturbation changes the morphology in a 
smooth manner without any nucleation of steps. There- 
fore it is possible to estimate the properties of films with 
a non- integer number of monolayers by interprolating the 
results of films of integer numbers of monolayers. Thus, 
a continuum model should adequately describe thin film 
growth behaviour for slightly vicinal surfaces. 

We assume that surface diffusion is the dominant mass 
transport mechanism. Gradients in the chemical poten- 
tial produce a drift of surface atoms with an average ve- 
locity, v, given by the Nernst-Einstein relation 



D s dfj, 
k R T ds' 



(1) 



where D s is the surface diffusion coefficient, s is the arc 
length, T is the temperature, ks is Boltzmann's constant 
and n is the chemical potential at the surface; i.e. it is the 
increase in free energy when an atom is added to the solid 
surface at the point of interest. Taking the divergence of 
the surface current produced by the atom drift gives an 
expression for the surface movement Jl9[| 



dh D s r)£l d d^i 
~di = k B T ds' 



(2) 



where rj is the number of atoms per unit area on the solid 
surface and Q is the atomic volume. 

In the continuum approximation fi = 04r-> where F is 
the free energy of the solid and SF/5h is the functional 
derivative of F. The free energy is composed of elastic 
and surface terms: 



F el + / dx jy/l + (dh/dxf 
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where 7 is the surface tension and F e i is the elastic free 
energy including any elastic contributions to the surface 
tension. 

In general the elastic free energy can be written as 
F e i = Fei'' + SFel, where F^' is the elastic free energy 
of a zero strain reference state. The elastic free energy 
can be written in terms of the elastic free energy density, 
f v , as F e i — J dxdyf v . f v is expanded as a power series 
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in the strain: f v = f u 
where fv is the free energy density in the zero strain 
reference state, 07? is the stress in the reference state and 
Cijki are the elastic coefficients of the material. In linear 
elasticity theory, deformations are assumed to be small 
and so terms of third order and higher are neglected. The 
stress-strain relationship is given by cry = J^-, which 
under linear elasticity gives Hooke's law: 



(0) . 
<?ij +Cijkiem- 



(4) 



We now briefly describe the equations which need to be 
satisfied in order to completely determine the equilibrium 
stress and strain in an elastic body. An elastic solid must 
satisfy the equations of mechanical equilibrium at every 
point in its interior (see e.g. |M ): 



djCTi 



6=0, 



(5) 



where £j is an external force on the solid. The solid must 
also satisfy the equations of equilibrium at its surface, 



IT, 



(6) 



where n is the exterior normal, and T n is the external 
force acting on the unit area surface element with normal 
n. As the strains are not independent but are linked 
via the displacements of the elastic body, they must also 
satisfy the equations of compatibility: 
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dkdie i3 + didjCui — d 3 die lk + did k e.j 



(J) 



Eqs. (j^),(0) and ^ along with the stress-strain re- 
lationships (jj) give us a system of equations, which are 
sufficient for the complete determination of the equilib- 
rium stress and strain in an elastic body. 

For the system we study, external body forces (e.g. 
gravity) are neglected. Hence Eq. (||) becomes 



djOi 



= 



for y < h(x). 



(8) 



Our system has periodic boundary conditions in the x 
direction and is infinite in the negative y direction. We 
shall assume that the forces on the upper surface due 
to surface tension (as given by Marchenko and Parshin 
pl[) are negligible in comparison to the forces due to 
the mismatch stress. This assumption is fulfilled as long 
as r ^ ^ £ wnere R is the radius of curvature of the 
surface and M is the plane strain modulus. For typical 
values of 7, M and e, this condition is satisfied when 
R is larger than the lattice constant. As typical surface 
features have length scales of the order of lOOnm this 
assumption is valid. Hence the boundary conditions are 
given by 



dij-rij =0 at y = h(x) 
Oij — > when y — > — 00 



(9) 



We now return to our discussion on the determination 
of the elastic free energy of both the reference state and 
the perturbed state. For each value of x, our reference 
state corresponds locally to a flat film of thickness h(x) 
constrained to have the lateral lattice constant of the 
substrate; i.e., = J dx dy f v °\h(x) , y) , where 

fv (h(x),y) is the elastic free energy density of a flat 
film of thickness h(x) with the substrate lateral lattice 
constant. We calculate the correction to the elastic free 
energy of the perturbed state, 8F e i, using linear elasticity 
theory. 

When looking at the stability of a strained flat film of 
thickness C, the obvious first choice for a reference state 
is that of a flat film of thickness C constrained to have the 
lateral lattice constant of the substrate. For later calcu- 
lations we must fully define the reference state and hence 
need to know its stress (C, y) and free energy density 

fv°\c, y). One simple approach to calculate these quan- 
tities would be to use linear elasticity with the unstressed 
film as a reference state. In linear elasticity a flat film 
of any thickness constrained to have the substrate lateral 
lattice constant and free to move in the y direction is in 
equilibrium and has the elastic free energy density of an 
infinitely strained film. Hence such a calculation does not 
predict any C or y dependence in and fv except for 
a step function at the film-substrate interface. For exam- 
ple in the case of plane strain where the mismatch strain 
is uniaxial (i.e., e xz = e yz = e zz = 0,e xy = 0,e xx = e), 



linear elasticity gives crj^ — Me and /i -* = Aie 2 /2, 
where M is the plane strain modulus. Therefore, varia- 
tion of the elastic free energy and stress of a flat film with 
film height is a nonlinear phenomenon, and a model out- 
side of linear elasticity theory must be used to calculate 
them. As will be shown in Section. Ill, small variations 
in the reference free energy density with film thickness 
are crucial in predicting wetting layer thickness, and a 
reference free energy density which has no variation with 
film thickness will lead to thin films that have no wetting 
layer. 

The disadvantage of our choice of the reference state 
is that the dependence of h on x leads to lateral vari- 
ations of the reference state. As a result, the reference 
stress does not satisfy the condition of mechanical equi- 
librium. However, the needed corrections vanish in the 
limit a I A — + 0, where a is the length scale over which 
stress varies in the y-direction and A is the lateral length 
of typical surface structures. This is because in this limit 
there are no lateral variations in the reference stress. As 
typical experimental islands have A ~ lOOnm, and as a is 
of the order of the lattice constant (see below), the cor- 
rections to the reference stress are small and have been 
ignored. 

Though linear elasticity cannot be used to calcu- 
late properties associated with the reference state, it 
can still be used to find the correction to the elastic 
free energy of the perturbed state, SF e i. For conve- 
nience we work in terms of the reference elastic free 
energy per unit length in the x-direction, f^\h(x)) = 

J^ x ^ dyf v a \h(x),y), instead of the free energy per unit 
volume. According to linear elasticity theory SF e i = 

In terms of the 



stress tensor, we find 



F d =fdxfM 



Jdx 



h(x) 



dy ( TiSijklCij&kl 



1 c (0) (0) 



(10) 



where we have used the inverted Hooke's law = 

Sijki{&ki — Cjy )■ Sijki are the compliance coefficients 
of the material. 

Using Eqs. (|J) and ( |Io| ) we arrive at an expression for 
SF/Sh 
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c - a (°) (°) 



h(x) 



y=h(x) 



d /l 



(11) 



where k is surface curvature, 9 is the angle between 
the normal to the surface and the y direction, and 
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7(0) = 7(0) + d 2 -f/d8 2 is the surface stiffness. The inte- 
grand in Eq. ([ll]) is the change in the linear elastic free 
energy density as the surface profile changes infinitesi- 
mally. Following Sokolnikoff [Q we now show that this 
term is second order in Sh in the absence of external sur- 
face or body forces and so can be neglected. When the 
surface profile changes, the strain in the body changes 
from eij to + e^-, where is of order Sh. The elas- 
tic free energy density can be written as f v + Af v = 



2 Cijkl i^-ij 



■e'ij){e k i 



z' kl ), and the change in elastic free 



energy density is Af v = c i]k ieije' kl + ^jkl^ij^kV Using 
Hooke's law (j|), the definition of strain and Eq. (|J), we 
rewrite Af v as Af v = dj(aiju' i ) + ^iu' i + \ciM e 'ij e 'u- Tne 
total change in the elastic free energy is J Af v dxdy — 
J T^u^ds + J ^u'idxdy + J ^c^kie'^e^dxdy, where the 
first term on the r.h.s. is an integral over the film sur- 
face. To obtain this we used Eq. (^). In the absence of 
external surface or body forces the first two terms on the 
r.h.s of the above equation vanish, and we are left with 
the equation J Af v dxdy — f \cijkie'^e' kl dxdy. is of 
order Sh, and hence the last term in Eq. ( pi] ) can be 
ignored for infinitesimal changes to the solid surface and 
we have 



SF _ 
-5h=^- 



dh 



1 c JP)JP> 



'ij u kl 



(12) 



y=h(x) 



As the above equation gives 4£ at the solid surface, all 
variables in the equation are also given at the surface. In 
particular crj? (h, y = h) is taken as the stress at the sur- 
face of a flat solid of height h(x) and hence must vanish 
when h < 0, since then the film is absent. df^?\h)/dh is 
determined by calculating how the reference elastic free 
energy of the solid changes as monolayers are added to 
the solid surface. When h < 0, df^'/dh = as the 
substrate is completely relaxed. In principle, Eq. ( |l2] ) 
should also contain derivatives of 7 with respect to h. 
However, we believe that the variation of surface tension 
with h away from a step dependence is due to elastic ef- 
fects. Since we included all elastic contributions in the 
zero-strain elastic free energy, we modeled 7 as a step 
function, taking the value of the substrate surface ten- 
sion for h < and the film surface tension for h > 0. 
Thus all partial derivatives of 7 with respect to surface 
height vanish and were omitted from Eq. ([l2]). 

Equations (||) and (|l2] ) form a complete model of sur- 
face evolution. In order to solve this model, the chemical 
potential (given by Eq. (|l2|)) for a given surface must be 
found, and so the linear elastic free energy density at the 
solid surface, j l ™ — ^Sijki&ij&ki\ y — h f x s must be calcu- 
lated. For an isotropic solid under plane strain with zero 
force boundary conditions the above expression simpli- 



fies considerably to give, /„ 



lin 
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2M 



where M is the plane strain modulus. Hence we must de- 
termine the stress at the surface of the film. This is done 
by finding the stress which satisfies both the linear elas- 
ticity equations (the equations of compatibility (^|) and 
equilibrium^)) and the boundary conditions (||). 

For an isotropic solid under plane strain, finding the 
stress which satisfies the linear elasticity equations and 
the boundary conditions can be reduced to finding the 
stress function, W, which satisfies the boundary condi- 
tions ( p|) a nd the biharmonic equation (see e.g. Timo- 
shenkol22| or Mikhlin [B3|): 



A 2 W = A(AW) 



d 4 W d A W 



d 4 W 



dx 4 dx 2 dy 2 dy 4 



0, (13) 



with a x:L 



<)-\Y 



and 



dx 2 



dy 2 , "xy — dxdy "yy 

In order to model the early evolution of faceted is- 
lands, and to study the effect of an anisotropic form 
of surface tension on the wetting layer, we used the 
cusped form of surface tension given by Bonzel and 
Preuss p3, which shows faceting in a free crystal: 7(0) = 
70 [1 + /3|sin(7r0/(20 o ))|], where (3 w 0.05 and O is the 
angle of maximum 7. The value of 70 was taken as 1 
J/m 2 in the substrate and about 75% of that in the film 
(as is the case for Si/Ge). This ensures a wetting layer 
of at least one monolayer. We considered a crystal which 
facets at 0°,±45° and ±90° with O = tt/8. The cusp 
gives rise to 7 = 00. However, a slight miscut of the low- 
index surface along the z direction leads to a rounding of 
the cusp, which can be described by 



7(0) = 70 l + /3Wsin^( 



'20n 



Q-2 



(14) 



where, for example, G = 500 corresponds to a miscut 
angle, A0 w 0.1°. As mentioned earlier all the results 
mentioned in this paper relate to surfaces with a very 
small miscut angle in the z direction. 



III. LINEAR STABILITY ANALYSIS 

In this section we carry out a linear stability analysis of 
Eq. (^) against a sinusoidal perturbation of wavenumber 
k, similar to that carried out in Ref. p5| for an infinitely 
thick stressed film. We thus look for a height profile of 
the form, h(x, t) = C + d(t) sin/err, which solves Eq. (Q) 
to first order in 8. To calculate the linear elastic energy 
we find the solutions of (|l3|), which satisfy the boundary 
conditions (^). a x °J vanishes because the film is hydro- 



statically strained, and a 



(0) 



since in the reference 



\y=h(xY 



state the force on the surface vanishes. Hence the only 
non-zero component of the reference stress is a x ° x \h,y). 
Stress functions of the form 

W = cr x °J (h, y)y 2 /2 + {A + By)e ky sin(kx) (15) 
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satisfy the biharmonic equation. This gives stresses of 
the form 

a xx = k[2B + {A + By)k]e ky sm{kx) + a^{h, y) 

G yy = -k 2 (A + By)e ky sin{kx) 

a xy = -k[B + (A + By)k]e ky cos(kx). 

To first order in 6 the stresses that satisfy the bihar- 
monic equation and the boundary conditions are given 
by: 

a xx = -5ke- kc cr x J(C, C){2 + (y - C)k]e ky sin(kx) 



'yy 



= k 2 e- kc a x °J(C, C)(y - C)e ky sm{kx) 

= 6ker kc a x °](C, C)[l + (y - C)k]e ky cos(kx). 



At the surface these stresses take the form: 

a xx = -2Ska x °J(C, C)sm(kx) 

+a x °J (C, C) + S sin{kx)da x °J /dh\ h=c 



'mi 



= 



a xy — Ska x °J(C, C) cos(kx). 

Note that all the derivatives in the Taylor expansions 
used in this analysis are with respect to h, the refer- 
ence film thickness and not with respect to y, the depth 
within the reference film. This is because in calculat- 
ing the chemical potential we are interested in how the 
free energy of the film changes as material is added or 
removed from the film surface; i.e. how the free energy 
of the film changes as the film thickness changes and not 
how the free energy density of the film varies within the 
film. 

Using the above stresses in Eq. (|l2|), we obtain the 
expression 



5F_ 
dh 



7« 



dh 

6 sin(fciz;) 

, 4f, (0) 



(0) 
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7^7 (cr^ + (Jyyf - ^— (a x °J + a. 



2M 

dh? 



2M 

2 



of . (^(C,C)) : , 2 
2k + 10 k 



dh 



(C, C), 



(16) 



where 70 = j(0 = 0). Combining the above equation 
with the evolution equation (|J) gives the following equa- 
tion for 5(t): 



d 4 = Kk 2 
dt 



where K 



-fc 2 7o + 2k 



D s rifl 2 
knT 



M 



Iff 
dh 2 



S, 



J h=C 



(17) 



Each term in the brackets in this 
equation has a simple physical significance. The first 



term is a surface tension term. Surface tension acts to 
reduce surface curvature, k, and so this term is negative, 
thereby reducing the perturbation amplitude, and is lin- 
ear in k ~ k 2 . The second term in this equation is a 
mismatch stress term. Regions of high stress have large 
chemical potential, and so atoms tend to detach from 
these regions and attach to regions of small chemical po- 
tential. In a mismatch stressed solid, valleys or cusps 
are regions of high stress, hence material moves from the 
valleys to the hills of a perturbed surface increasing per- 
turbation amplitude. The contribution of this term is 
propotional to the density of valleys, which is linear in k. 
The last term is a reference state term. If d 2 /dh 2 > 
it costs more energy to add a monolayer to a flat film 
than to remove a monolayer, and hence it costs energy 
to perturb a film. Thus, positive d 2 ff/dh 2 stabilizes a 
flat thin film, whereas negative d 2 ff/dh 2 leads to an 
instability. Obviously this reference state term is present 
even if the film is flat and hence is independent of k. 

Equation ( |l7|) implies that the flat film is stable at all 
perturbation wavelengths as long as 



a x °](C, C) 



M 2 



< 7o 



d 2 f, 



(0) 
el 



dh 2 



(18) 



h=C 



and the equality holds at the critical wetting layer 
thickness, where perturbations of wavenumber k — 

fn) ~\ 2 

g xx (C, C) j{M 70) are marginal. 70 is positive if 9 = 

is a surface seen in the equilibrium free crystal ^6|. As 
mentioned earlier 70 — > 00 at a perfect facet, and is large 
and positive on a surface with a small miscut, as is the 
case for most of the materials used in epitaxial films. 
Therefore, a linearly stable wetting layer of finite thick- 
ness can exist only if d 2 jfjdh 2 > 0. Note that for 
the wetting layer to have a finite rather than an infi- 
nite thickness, d 2 /dh 2 must decrease to a value less 
than the l.h.s. of Eq. (|lj) as the thickness of the film 

increases. a x °J (C, C) depends linearly on the lattice mis- 



match e, and hence the l.h.s. of (18) is proportional to 
e A , while the r.h.s. of ( |l8| ) is proportional to e 2 due to 
the dependence of on lattice mismatch. Therefore, 
if d 2 ff/dh 2 > 0, the thickness of the wetting layer in- 
creases with decreasing lattice mismatch and diverges in 
the limit e — > 0. 

Note that the maximum thickness of a flat film which is 
stable to infinitesimal perturbations is given by ( [T^ ) when 
the equality holds. A film slightly thicker is unstable to 



a {0) 

u XX 



(C,C) 



perturbations of wavelength A = 2nMjo I 
For films which are nearly perfect facets with small mis- 
cut angles these wavelengths are larger than the typi- 
cal sample size and so practically such perturbations will 
never occur. However as will be explained in Section 
V, the film can be nonlinearly unstable to smaller wave- 
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length perturbations of a non-zero amplitude at phys- 
ically reasonable wavelengths. Hence the inequality in 
( |l8| ) is only useful in predicting the stability of films with 
large miscut angles or above the roughening temperature. 
At small miscut angles the stability of the film to large 
perturbations will predict its maximum thickness. This 
issue is discussed in more detail in Section V. 

IV. CALCULATION OF THE NONLINEAR 
ELASTIC FREE ENERGY OF A FLAT FILM 

As can be seen from both Eq. (|l^ ) and Eq. (|l8|). the 
dependence of the nonlinear elastic free energy of a flat 
film feP(h) on film thickness, h, is vital in order to de- 
termine both wetting layer thickness and thin film evolu- 
tion. This free energy depends strongly on the mismatch 
stress Oy , and its dependence on the y coordinate. As 
a result of the sharp interface between the substrate and 
the film, we expect to behave as a step function of 
y with small corrections due to elastic relaxation. If wc 
ignore these small corrections, the resulting free energy 
f^\h) is proportional to film thickness, and its second 
derivative vanishes. Hence according to Eq. (|l8|), the 
thickness of the critical wetting layer vanishes. The cor- 
rection due to elastic relaxation is therefore extremely 
important. As discussed earlier, this correction vanishes 
within linear elasticity theory. This led some investiga- 
tors [^7j] to claim that the variation in free energy over the 
interface was due to nonelastic effects, e.g. film-substrate 
material mixing over the interface. However, we claim 
that this is not necessary, since nonlinear elasticity can 
explain the corrections to the step-function form of the 
free energy. 

The nonlinear clastic free energy of the reference state, 
fe?\h), is calculated for a solid with a flat surface of 
height h. Hence in order to calculate df^ /dh, we 
determine for flat solids of heights h + 5/2 and 
h-5/2, and use the estimate df^/dh = [f^\h + 5/2) - 

fj?P(h — 5/2)]/ 5. Ideally, first principles, substance spe- 
cific calculations should be performed in order to evaluate 
&xx {h, h) and (h) , and we intend to carry out such 
calculations. However, the qualitative general behaviour 
of (h) can be obtained from much simpler models. To 
demonstrate this point we carried out the calculation for 
two-dimensional networks of balls and springs of varying 
lattice- type and spring constants. In these calculations 5 
is one monolayer, and (h) is calculated at film thick- 
nesses of integer numbers of monolayers from no film up 
to 10 monolayers of film. Values of (h) for film heights 
of fractional monolayers are interprolated from the values 
calculated at integer monolayer heights. 

In the ball-and-spring model the balls are connected 
by springs that obey Hooke's law. Note that this does 



not imply that the stress-strain relationship of the ball- 
and-spring network is linear (a discussion of this point 
can be found in Feynman's Lectures p3]). The natural 
spring length has a step variation over the interface, and 
the balls are placed on a lattice with the substrate lat- 
tice constant. Thus balls in the substrate are connected 
by springs of their natural length, whereas balls in the 
film are connected by springs that have undergone a hy- 
drostatic transformation strain and have length larger 
than their natural length by a factor of 1 + e, where e 
is the homogeneous strain in the film. The network was 
then allowed to relax, with the film free to move in the 
y-direction, and periodic boundary conditions being ap- 
plied in the cc-direction to ensure that the system bound- 
aries in this direction were fixed to the natural substrate 
length. 

We calculated the mismatch stress a x °x within the re- 
laxed film and at the film surface. We also calculated 
the dependence of mismatch surface stress a x °J (h, h) and 
the nonlinear elastic free energy /S° (h) on varying film 
thickness, h. We carried out these calculations for various 
two dimensional networks of balls and springs with vary- 
ing spring constants. An example of the ball-and-spring 
model on a fee- like lattice is shown in Fig. (0). 



y 




one 
monolayer 



x 

FIG. 2. Example of a fcc-like lattice. The circles represent 
balls. The curvy lines are springs of spring constant ki, and 
the dashed lines represent springs of spring constant k%. 

Simulations showed that whilst individual atoms were 
free to move in the x-direction, they actually moved only 
in the y-direction. The relaxation in the y-direction de- 
pended on the film thickness and on the depth of the 
atom in the lattice but was independent of a;. In general 
balls at depth of more than 3 monolayers into the sub- 
strate experienced no stress. The stress experienced by 
balls close to the interface depended on the lattice type, 
spring constants and ball position within the monolayer. 
A few monolayers into the film, balls experienced the 
stress of an infinite thickness film, Me. At the film sur- 
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face balls showed large relaxation. Figure || shows an 
example of the mismatch stress, &xx (h, y) in a typical 
fcc-like lattice. 
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FIG. 3. Mismatch induced stress, ai°J(h,y)/Me, in a typ- 
ical fcc-like lattice. The film is 8 monolayers thick and the 
substrate is 10 monolayers thick. h m i is the thickness of one 
monolayer. 

Note that the springs in a simple square lattice could 
relax completely in the y-direction. Therefore, for such 
a lattice the relaxation is independent of spring depth 
within the film or film thickness, and (h) varies lin- 
early with film thickness. Hence for a square lattice 
d 2 fei\h) I dh 2 — 0. Only when diagonal bonds, such as 
those in a fee lattice were present did the springs show 
depth dependent relaxation such as that described above. 
The inability of the springs to completely relax due to the 
presence of diagonal bonds was a necessary condition for 
fei^Q 1 ) t° vary nonlinearly with film height. In such in- 
completely relaxed films the nonlinear dependence of ffy 
on h arises from the elastic relaxation at the surface and 
its coupling to the relaxation at the film-substrate inter- 
face. A similar effect should occur in real systems due to 
surface reconstruction, for example. 

A typical behavior of df^ /dh is shown in Fig. |4|, where 
it is seen that f£ (h) indeed depends on the thickness h. 
Moreover, the model predicts that d 2 f^/dh 2 > and 
decreases with increasing film thickness, and therefore 
according to the inequality (18) and the discussion fol- 
lowing it, there should be a linearly stable wetting layer, 
whose thickness is finite and increases with decreasing 
lattice mismatch. 



Whilst the detailed dependence of df^\h)/dh on 
film thickness close to the substrate-film interface [< 
3 monolayers] varied between different networks, it 
showed the same general behaviour. In all systems 
d 2 fei (h)/dh 2 showed exponential decay with a decay 



length of about a monolayer from the interface. The 

was independent of lat- 



dimensionless quantity ^ 



dh 



tice mismatch sign or magnitude, df^ j dh increased with 

film thickness. As the film thickness increases df^ /dh 
asymptotically approaches the elastic free energy density 
of an infinite film, Me 2 /2, as expected. 
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FIG. 4. Variation with film thickness of the elastic free en- 
ergy of a relaxed ball-and-spring system, df^ /dh, as a func- 
tion of film thickness h. The free energy is normalized to the 
infinite film linear elastic energy density, ^Me 2 . h m i is the 
thickness of one monolayer. 

We can explain the increase of df^ / dh with film thick- 
ness intuitively. df^P / dh is the change in the free energy 
of a film as a monolayer is added. When a monolayer 
is added to a thick film the contribution to the free en- 
ergy from the surface atoms and the interface remains 
the same and the energy added is effectively that of a 
monolayer in the 'bulk' of the film. Atoms in the bulk 
of a thick film are more constrained than those in a thin 
film where the atoms are relatively free to move and re- 
lax. When a monolayer is added to a thin film the energy 
change is in between that of the constrained thick film 
bulk atoms and the relaxed surface atoms. Hence more 
free energy is needed to add a monolayer to a thick film 
and df^/dh increases with film thickness. 

For the calculations used later in this paper we used 
the function 



df^(h)/dh = ^£-[1 - 0.05exp(-tyftr. 



for h > 0, 
(19) 



and dfel (h)/dh = for h < 0. h m i is the thick- 
ness of one monolayer. The exponential decay form of 
df^\h)/dh fits the results given by Tersoff p3| . In pre- 
vious works |ll],[l2|,^7j on the physics of the wetting layer 
it was assumed that the reference state energy variation 
is a smooth function of h, mainly in order to avoid non- 
analyticities at the interface. In contrast, our reference 
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state energy variation behaves as a step function of the 
surface height with a small but important correction. 

The deviations in the mismatch stress (averaged over 
the surface monolayer) at the film surface cr^x (h, h) from 
a step function, a x °x (h, h) = Me when h > 0, was shown 
to be small (< 5%) but dependent on spring constants 
and lattice type (See Fig. (||)). As variations in a x °J (h, h) 
only slightly alter the wetting layer thickness predicted 
from Eq. jig), we decided to use the step function form 
of mismatch stress. 
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(b) 
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FIG. 5. Variation of mismatch induced stress at the film 
surface, a x °J (h, h) /Me, for varying film thickness, h. The 
lattice is fee- like (see Fig. |^). The dashed lines with circles 
denote a x °J /Me of atoms at the top of a monolayer. The solid 
lines with circles denote a^x /Me of atoms at the bottom of 
a monolayer. The solid lines without circles represent the 
averages over a monolayer. The spring constants 

which were used are (a) fei = 1 and ki = 1 and (b) ki = 1 
and k'2 = 5. 

Combining the behavior of df^ / dh from Eq. ([l9]) with 



the inequality (|18|), we obtained an expression for the 
linear stability wetting layer thickness, h c : 



h c /h ml = max {l, In [y /(AOM£ 2 h m i)] } 



(20) 



Thus, the wetting layer thickness increases with decreas- 
ing lattice mismatch, as observed in experiments. 



V. NUMERICAL SIMULATIONS OF THIN FILM 
GROWTH 

As explained at the end of Section III, for systems with 
small mismatch and/or small vicinality (miscut angle) 
one has to go beyond linear stability in order to under- 
stand their evolution. To this end we carried out nu- 
merical simulations of the evolution of the strained film. 
The evolution equation (^) given by Mullins [I9|, which 
is derived from the Nernst- Einstein relation |l ) , includes 
derivatives of the chemical potential, n, along the surface. 
This chemical potential is defined as the change in free 
energy when an atom is added to the surface. Continuum 
theory assumes that the free energy change when the sur- 
face is changed by an infinitesimal amount is proportional 
to this chemical potential (jj, = f^f?)- However, when we 



solved the evolution equation by directly calculating the 
chemical potential from Eq. ( |l2| ) at points along the film 
surface we experienced numerical instabilities. 

We have come up with the following solution to this 
problem. The Nernst-Einstein equation can be derived 
by considering material of atomic volume moving along 
the solid surface. When material jumps between neigh- 
boring atomic sites, it must cross a free energy barrier 
of AF± = E d + (F± - F )/2, where E d is the potential 
barrier for diffusion, F± is the free energy of the film af- 
ter material has been moved, and Fq is the free energy 
of the film before material is moved. The positive and 
negative signs stand for forward and backward jumps, 
respectively. This leads to the following equation for ma- 
terial velocity along the surface 



E d /k B T ( -(F + -F )/2k B T _ -(F--F„)/2k B T 



D, 



(e 

±^ e -{F+~F )/2k B T _ e -(F--F )/2k B Tj 



(21) 



where u is the attempt rate and a is the jump length. 
When AF± is small, this equation gives the Nernst- 
Einstein relation (Q). 

We solved Eq. (|2l]) using the following numerical 
scheme. For every two adjacent points on the surface, the 
surface height of the left point was changed by ±Sh and 
of the right point by =p<5/i so as to give a transfer of mate- 
rial of atomic volume backwards and forwards along the 
surface respectively. The change in surface free energy, 
J dx 7yl + (dh/dx) 2 , was calculated for this material 
transfer. The change in the elastic free energy was calcu- 
lated at each point using the integral of Eq. (fflh, 5F = 



± 



dh 



^Sijki&ijCrkl 



2^ijkl <J ij 



(0) to) 

°kl 



8h. 



y=h(x)_ 

The linear elastic energy was calculated by solving the 
biharmonic equation (|l^) with the boundary conditions 
(^) . This was done by solving a boundary integral equa- 
tion in terms of the complex Goursat function, the details 
of which can be found in the paper of Spencer and Meiron 

M. 



A. The stability of thin Alms 

According to Eq. (^p|), anisotropic surface tension 
greatly enlarges the linearly stable wetting layer thick- 
ness. Does this conclusion survive beyond linear stabil- 
ity analysis? When a linearly stable flat film is perturbed 
strongly so that the surface orientation in some regions 
is far from the 8 — direction, the local surface stiffness 
in these regions is much smaller than the 8 = stiff- 
ness. This tends to destabilize the linearly stable film. 
Indeed, we carried out numerical simulations (using the 
procedure described above) that showed that films thin- 
ner than the linear wetting layer were unstable to per- 
turbations greater than a certain critical amplitude (see 
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Fig. ^) . Hence films thinner than the linear wetting layer 
thickness are metastable. When large perturbations were 
applied, faceted islands developed in the film, which un- 
derwent ripening at later stages of the evolution. 
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FIG. 6. Evolution of a randomly perturbed film, in which 
perturbations were larger than the critical perturbation am- 
plitude. Lattice mismatch in this film is 4%. The initial film 
surface is shown as a thin solid line. The dashed line shows the 
film surface at a later time. The linear wetting layer thickness 
is shown as a thick solid line. 

We carried out simulations for films perturbed by ran- 
dom perturbations and by perturbations of a single wave- 
length. The critical perturbation amplitude, 8 C , depends 
on the wavelength of the perturbation, A, taking its min- 
imal value 



min <5 C (A) 



(22) 



at X/l Q ~ 10 — 50, where l a = 2j Q /Me 2 . 5™ in monolayers 
is plotted as a function of lattice mismatch in Fig. [?]. The 
linear wetting layer thickness for G = 500, M = 1.5 x 
lQ ll N/m 2 and h m i — 5A is also shown for comparison. 

<$™ was found to be proportional to e~ 2 . The e~ 2 de- 
pendence is expected for an infinite film as in this case the 
evolution equations (Eq. (||) together with Eq. ([12])) can 
be made spatially dimensionless by scaling all lengths by 
Iq. Hence all perturbations of size 8/Iq and with the same 
dimensionless wavenumber Hq will evolve identically. 

<5™ was largely independent of cusp smoothness G, un- 
like the linear wetting layer thickness which depended 
strongly on G. This suggests that unlike the linear wet- 
ting layer thickness the critical perturbation amplitude 
can be used in predicting the outcome of experimental 
thin film growth. The mean square amplitude of the 
random perturbation needed to destabilize thin films was 
also shown to be largely independent of G and was pro- 
portional to £~ 2 . 
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FIG. 7. Variation of the minimal critical perturbation am- 
plitude, <5™, and linear wetting layer thickness, h c , with lattice 
mismatch, e. The minimal critical perturbation amplitude, 
8™/hmi is represented by the solid line. The linear wetting 
layer thickness, h c /h m i, is represented by the dashed line. 
The dotted line shows the size of one monolayer, h m l, for 
comparison. 

When the lattice mismatch is small, 8™ is much larger 
than the linear wetting layer thickness (see Fig. [?]). 
Hence, flat films thinner than the linear critical thick- 
ness are stable at small lattice mismatch. As the linear 
critical thickness at small lattice mismatch is very large, 
we expect the film to first become unstable to misfit dis- 
locations. This in indeed seen in experiments J^J^]. 

At intermediate lattice mismatch, <5™ is of the order of 
a few monolayers. Hence we expect such a film to become 
unstable as perturbations of this amplitude are physically 
likely. In this regime films should develop growing per- 
turbations at wavelengths given by X/l a ~ 10 — 50. This 
corresponds to wavelengths of a few hundred nanometers. 
This typical wavelength decreases as lattice mismatch in- 
creases, agreeing with experiment ||||. As e increases, 
8™ decreases in this regime from about 10 monolayers 
to approximately one monolayer, we expect the thick- 
ness of the film needed to support such perturbations to 
correspondingly decrease. Such a trend is seen in exper- 
iments 1^^] • In order to compare quantitative wetting 
layer thickness and its lattice mismatch dependence with 
experiments we will have to carry out first principles, 
substance-specific calculations to evaluate f^\h). How- 
ever, the general qualitative trends predicted here agree 
with experimental observations. 

Looking at Fig. [?] we see that at intermediate lattice 
mismatch the critical perturbation amplitude, 5™, and 
the linear wetting layer thickness are of the same order 
of magnitude (several monolayers). This could mean that 
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we have to also consider the linear wetting layer thickness 
with its strong dependence on the surface miscut angle 
when deciding whether or not a thin film will be un- 
stable. However, the infinitesimal perturbations needed 
to perturb the linear wetting layer occur at wavelength, 
A = 



where we have assumed the reference state mismatch 
stress is given by a step function, Uxx {h, h) = Me when 
h > 0, and a^(h,h) — when h < 0. Using the gen- 
eral form (obtained from the ball-and-spring model) of 



2n 
Ms- 



d 2 f^/dh 2 



A2L 



exp(—h/h m i), where x is a constant, 



7o , whereas the typical wavelength at which crit- giveg the f li ow i n |' so l u tion for perturbation growth: 



ical amplitude perturbations first appear in the film is, 
A ~ 10j^j7o- Using expression (|l4|), the ratio between 
these two wavelengths is approximately equal to G. For 
physical values of G the linear wetting layer perturba- 
tions will have wavelengths of the order of 100/im which 
is larger than the typical sample size, whereas the wave- 
length which corresponds to the minimal critical pertur- 
bation is much smaller (~ lOOrtm). Therefore physical 
thin films should first become unstable when the film 
thickness is large enough to support perturbations larger 
than the critical perturbation amplitude. 

For very large mismatch, a perturbation smaller than a 
monolayer is sufficient in order to destabilize the linearly 
stable wetting layer. Therefore, in practice, the wetting 
layer is a single monolayer in this case. 



VI. EARLY EVOLUTION OF THIN FILMS WITH 
MATERIAL DEPOSITION 

We carried out our calculations with two different 
types of material deposition: The first type is deposition 
at a steady rate in the vertical y-direction, corresponding 
to any directed deposition (eg, molecular beam epitaxy). 
The evolution equation (^) then becomes 



dh D s r]£l d d[i 
~di ~ k B T dx'ds 



Vd, 



(23) 



where Vjj is the material deposition rate. 

The second type is deposition constant in the direction 
perpendicular to the film surface, corresponding to liq- 
uid phase epitaxy, for example. Early growth with this 
method of deposition has been studied by Chiu and Gao 
PJ . In this case the evolution equation becomes 



dh D s r]D, d d^i Vd 
dt ksT dx ds n y 



(24) 



where n y is the y-component of the normal vector to the 
surface. 

We performed linear stability analysis in order to ob- 
tain the analytical early evolution equation of a per- 
turbed thin film. This analysis is valid for both types of 
material deposition. Under steady deposition, Eq. (|lT 
becomes 



dS_ 

dt 



K 



-k% + 2k 3 Me 2 - fc : 



d 2 fS\c + V D t) 
dh 2 



5(t) = S exp {K [{-k% + 2k 3 Me 2 )t+ 
k 2 X (Me 2 /2V D )exp(-C/h ml )(exp(-V D t/h ml ) 



1)]}. (26) 



Note that in linear stability analysis a perturbation in 
an infinite film decays when k > 2Me 2 /jQ and grows 
exponentially when k < 2Ms 2 /jq. 

For isotropic surface tension, numerical computations 
showed that when k* < k < 2Me 2 /"jo, with k* w 
0.875 x 2Me 2 /jo, both methods of deposition lead to 
cusp formation in the surface valleys. The cusps initially 
evolve according to the linear evolution equation ( p6|) and 
then slow and reach a steady state morphology. Spencer 
and Meiron |2£j observed such steady states in infinitely 
thick stressed films. However when k < k* , surface evo- 
lution depends on the method of material deposition. 
When deposition is constant in the vertical y-direction 
increasingly sharp cusps form in the surface valleys (see 
Fig. ||a), which continue to grow exponentially. In con- 
trast when deposition is constant perpendicular to the 
surface at very high deposition rates cusp formation is 
slowed (see Fig. ||b) and the surface shows signs of reach- 
ing a steady-state morphology as for k > k* . 
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FIG. 8. Film evolution at very high deposition rates 
(Vd = 1200nm/s) when surface tension is isotropic. 
k = 0.75 x 2A/e 2 /7o. (a): Deposition is constant in the verti- 
cal j/-direction. (b): Deposition is constant perpendicular to 
the film surface. 

This can be seen in Fig. ^. The plot (shown as squares 
in the figure) starts as a graph of constant positive slope 
representing an exponentially growing perturbation as 
predicted by linear analysis. However this growth slows 
and the graph approaches the flat line representative of a 
steady-state morphology. Note that in comparison when 
deposition is constant in the vertical y-direction (shown 
as circles in Fig. [)]), the film evolves according to the lin- 
ear evolution equation. Under deposition perpendicular 
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to the film surface, when a cusp begins to form, mate- 
rial is deposited more rapidly on the steep cusp slope, 
hence slowing cusp formation. When deposition is con- 
stant vertically it only effects surface evolution indirectly 
by raising the average surface height. Though deposi- 
tion rates of this magnitude are not generally used in 
experiment it is nevertheless physically interesting to ob- 
serve the difference in surface evolution between the two 
growth methods. 
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FIG. 9. Film evolution at very high deposition rates 
(Vd = 1200nm/s) when surface tension is isotropic. The 
line represents the linear evolution equation, the circles rep- 
resent the results from the numerical simulation when depo- 
sition is constant in the {/-direction and the squares represent 
the results from the numerical simulation when deposition is 
constant perpendicular to the surface, k = 0.75 x 2Me 2 /70. 
The parameters used here are: D a — 3.599 x 10~ 13 m 2 /s, 
Q = 1.38 x l(T 29 m 3 , r] = 1.74 x 10 19 m~ 2 , T = 700K, 
k = 10 8 m-\ 70 = 70 = 1, e = 2%, M = 1.67 x 10 11 , x = 1, 
C — 0.75 monolayers and h m i = 2nm. Note the deviation 
from the linear stability analysis results when deposition is 
perpendicular to the surface. 

When the deposition is constant in the vertical y- 
direction, the film evolves according to the linear evo- 
lution equation, even after the surface is no longer a sine 
function and cusp formation occurs. This can be seen in 
Fig. |l^ which compares results from the numerical sim- 
ulation with the results predicted by the linear evolution 
equation (p6|). Figure [| shows a very clear cusp forma- 
tion in the surface morphology, whilst for the same time 



Fig. [n] shows the sharp cusp growing only slightly faster 
than predicted by linear analysis. This slight deviation 
is expected as the stress in a cusp valley is larger than in 
a sine valley hence accelerating perturbation growth. 

When the surface tension is anisotropic the surface evo- 
lution is very different from that predicted by the linear 
analysis. As can be seen in Fig. [l^, a perturbation in 
an isotropic film decays until the film surface reaches a 
height at which the film is linearly unstable to pertur- 
bations at that wavelength . No matter how large the 
deposition rate, at any given time a perturbation in a 
thin film is always smaller than a perturbation of the 
same initial size in an infinite film due to the finite time 
spent in the linear wetting layer. 
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FIG. 10. Film growth when surface tension is isotropic and 
deposition is constant in the {/-direction. The symbols repre- 
sent the numerical simulation results and the lines the linear 
evolution equation values, k = 0.75 x 2Me 2 /^q. The param- 
eters used are the same as those in Fig. |^. 

When surface tension is anisotropic, initially the per- 
turbation amplitude decreases as the surface facets. The 
rate and amplitude of the decrease is independent of ei- 
ther C or Vd- The film then grows or decays depending 
on whether the perturbation amplitude S(t), is larger or 
smaller than the critical perturbation amplitude S c (k, h), 
at h = V D t + C. When S > S c (k, h = C) the per- 
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turbation grows immediately after faceting. When So < 
S c (k, h = C) the perturbation initially decays though it 
will start to grow again if the deposition brings the film 
surface to a height at which 8{t) > S c (k, h — Vot + C). 

We now turn to characterize the evolution of the film 
in terms of all the relevant physical variables. There are 
five independent variables that can effect the evolution: 
Vd, C, k, t and e. In addition, there are four relevant con- 
stant parameters: h m i,K,j and M. We can replace the 
five independent variables by the dimcnsionless variables: 
V D t/h mh C/h mh Kj kH, KMe 2 k% KMe 2 k 2 t/h ml . 
The idea is that the description of the evolution of a film 
with isotropic surface tension becomes simpler in terms 
of these variables. This becomes clear when we look at 
how a perturbation in a thin film grows in relation to 
the growth of a perturbation of the same initial size in 
an infinite film. Quantitatively, this is described by the 
relative perturbation height, Sr — S(t,C)/S(t,C = oo). 
As can be seen from Eq. (p6|), Sr depends only on three 
of the five independent dimcnsionless variables: 



Sr = e 



Kk- 



-exp(-C/h mt )(ex P (-V D t/h ml )-l) 



(27) 



This can be summarized in the scaling law: 



5 R (V D ,C,k,t,s) = 5 R (V D t/h mll C/h ml ,KMe 2 k 2 t/h ml ), 

(28) 

which implies that Sr depends only on three of the five 
scaling variables. The manifestation of this scaling be- 
haviour is data collapse. For example, when Vd = and 
C is fixed, all the curves of <5^(i) for different values of k 
and e fall onto a single curve if plotted as a function of 
k 2 e 2 t rather than t. 

Does this scaling law survive beyond linear analysis? 
We looked for scaling when deposition was constant in 
the vertical y-direction for both isotropic and anisotropic 
surface tension. As mentioned earlier when surface ten- 
sion is isotropic the film continued to evolve according 
to the linear evolution equation (^6|) long after it left 
the linear regime and hence the scaling relation (^8|) also 
held. 

Growth under anisotropic surface tension is very dif- 
ferent from that given by the linear evolution equa- 
tion (^6|), and hence the scaling relation (28) does 
not hold. Does this mean that the physics of 
anisotropic surfaces is more complicated and depends 
on all five independent variables? It turns out the 
answer to this question is no to a good approxima- 
tion. To see this we define five generalized dimcnsion- 
less variables: VDt/h m i, C/h m i, K^$k A t, KMe 2 k 3 t, and 
K7^/*+i M -P-q/2 e -2p- qm / h i/*+4 w when p = 2 an d 

q = — 6 we regain the dimensionless variables govern- 
ing the linear evolution equation. We found numerically 
that in the case of anisotropic surface tension, Sr approx- 
imately obeys the scaling law: : 



S R (k,V D , e,t,C) 



SR{Vr>t/h m i, C/h 

Till • 



p+g/2+l n, r -p-q/2 £ -2p-q ^ //j<j/ 2 + 4 



) (29) 



with p ~ 2.37 and q ~ —6.5. Again, Sr depends only 
on three of the five scaling variables, which implies data 
collapse. This relation was very robust. We verified it 
for different G, variation of k of an order of magnitude, 
variation of e by 100% and deposition rates of between 
and 120000A/s. Figures [II] and [l2] show this scaling in 
the form of data collapse when C = 2ML. Data collapse 
when Sr is plotted against a variable proportional to the 
third scaling variable in Eq. (|2^) can be seen in Fig. [ll]. 
Here there is no deposition, k varies by over an order 
of magnitude and e by 100%. As can be seen the data 
collapse is not exact but holds to a good approximation. 
Figure [l2| shows data collapse when Sr, is plotted against 
variables proportional to the first and third scaling vari- 
ables in Eq. (p9|) . Deposition rates vary by six orders of 
magnitude, k varies by over an order of magnitude and 
e by 100%. Data collapse is shown by all curves falling 
onto a single surface. 
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FIG. 11. Scaling of relative perturbation height for dif- 



ferent k (k = 8 x 1CT 5 ^ -» 2 x lO -3 ^-) and e 
(e = 4% — » 7.5%) with zero deposition. p = 2.37 
and q = —6.5. Note that the variable used to plot this 
graph 2 ~ m ° / 04 ' i k p e~ 2p ~ q t is proportional to the scaling vari- 
able Kj p+q/2+1 M- p -^ 2 e- 2p - 1 k p t/h q ^ +i . 

The scaling relationship (^), however, only held when 
the initial perturbation 5q was larger than the critical 
perturbation at that wavenumber, k, and initial film 
height, C; i.e for So > S c (k,C). This is probably be- 
cause when So > S c (k, C) perturbations of a thin film 
and of an infinite film evolve similarly. Both perturba- 
tions initially decay whilst faceting and then continue to 
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grow. On the other hand, when S c (k, oo) < So < S c (k, C) 
a perturbation of a thin film decays whereas an infinite 
film perturbation facets and grows. In this regime scaling 
laws were not found. When 5q < S c (k, oo) the perturba- 
tion decays in both the thin and infinite film. When 5q 
is small enough the scaling laws (|2|) derived from the 
linear evolution equation are regained. 




FIG. 12. Scaling of relative perturbation height for differ- 
ent k (k = 8 x 1(T 5 ^- -> 2 x 1(T 3 ^-) and e (e = 4% -» 7.5%) 

V 70 70 J v y 

and V D (Vd = 0, 0.12, 1.2, 12, 120, 1200, 12000, 120000A/s). 
p = 2.37 and q — —6.5. Note that the variables used to 
plot this graph 2 ^j^-k p e~ 2p ~ q t and Vbt are proportional 
to the scaling variables Ky p+q/2+1 M- p - q/2 £ - 2p - q k p t/h 9 ^ +4 
and V D t/h m i. 
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